### This script will run class-stat on ASCII files and summarize the outputs
### Also, it will produce a bunch of plots comparing the stats between models within species, for all species

# Load library

library('SDMTools')

## Identify directory for ASCII files

asc.dir = '/home1/99/jc152199/MAXENT/output/'

# Identify species of interest

species = c('SAPCZEC','SAPTETR','CARRUBR','LAMCOGG','GNYQUEE','SAPBASI','LAMROBE')

## Blank object for data binding within loop

outdata = NULL

### Loop through species, extracting and summarizing ClassStat data

for (spp in species)

	{
	
	### Loop through different models within species, extracting and summarizing data
	
	for (mod in c('BIOCLIM','microCLIM'))
	
		{
		
		### Read in the ASCII
		
		t.asc = read.asc.gz(paste(asc.dir,spp,'/',mod,'/output/',spp,'_',mod,'_Realized.asc.gz',sep=''))
		
		### Calculate total summed ES
		
		ES=sum(t.asc[which(is.finite(t.asc)==T)])
		
		### Apply threshold to create a binary matrix
		### Modified for proportional ES values
		
		t.asc[which(t.asc>0)]=1
		
		### Run class stat

		t.stat = ClassStat(t.asc,cellsize=250,bkgd=NA,latlon=TRUE)
		
		### Extract data of interest
			
		t.out = data.frame(spp=spp,mod=mod,t.stat[2,],ESsum=ES)
		
		### Bind to a dataframe
		
		outdata = rbind(outdata,t.out)
		
		### Report progress
		
		cat('\n',spp,' ',mod,' Summary Completed\n',sep='')
			
		}
		
	}
	
# Done

### Write outdata

write.csv(outdata,file=paste('/home1/99/jc152199/MAXENT/ClassStat_RealizedDist.csv',sep=''),row.names=F)
outdata = read.csv('/home1/99/jc152199/MAXENT/ClassStat_RealizedDist.csv',header=T)

### Establish a directory to write plots into


plot.dir = '/home1/99/jc152199/MAXENT/images/'

### Perform I-statistic similarity calculations within species between models and compare raw frequency distributions within species between models

### Istat data

istatdata = NULL

for (spp in species)

	{
	
	### Read in the ASCIIs
		
	bcmap = read.asc.gz(paste(asc.dir,spp,'/BIOCLIM/output/',spp,'_BIOCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	acmap = read.asc.gz(paste(asc.dir,spp,'/microCLIM/output/',spp,'_microCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	
	### Density
	
	bc.density = density(bcmap[which(is.finite(bcmap)==T)])
	ac.density = density(acmap[which(is.finite(acmap)==T)])
	
	### Plot density
	
	png(paste(plot.dir,spp,' - Density.png',sep=''),height=18,width=18, units='cm',res=1000)
	
	### Plot
	
	plot(bc.density, xlab='ES', ylab = 'Freq.',main=spp, col='red')
	
	### Lines
	
	lines(ac.density,col='blue')
		
	### Caluclate I-statistic
		
	istat = Istat(bcmap,acmap)
	
	### Summarize species data
	
	tout = data.frame(spp=spp,Istat=istat)
	
	### Bind to blank object
	
	istatdata = rbind(istatdata,tout)
	
	#### Put on a legend
	
	legend('topright', legend=c(paste('IStat - ',istat,sep='')), text.col='black')
	
	#### Shut device driver
	
	dev.off()
	
	#### Report progress
	
	cat('\n',spp,' - Completed\n',sep='')
	
	}
	
# Done

############### Plotting difference between accuCLIM and BIOCLIM MaxEnt distributions

for (spp in species)

	{
	
	### Read in the ASCIIs
		
	bcmap = read.asc.gz(paste(asc.dir,spp,'/BIOCLIM/output/',spp,'_BIOCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	acmap = read.asc.gz(paste(asc.dir,spp,'/microCLIM/output/',spp,'_microCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	
	### Calculate the difference
	
	dmap = acmap-bcmap
	
	### Convert ASCII to vector
	
	cmap = c(dmap)
	
	### Range of values in difference ASCII
	
	lims = round(range(cmap,na.rm=T),2)

	### Create a custom color palette for plotting

	### Start by establishing break points for colors

	b1 = sort(-(seq(.01,abs(lims[1])+.01,.02))) # Breaks below zero
	b2 = seq(.01,lims[2]+.01,.02)	# Breaks above zero
	b3 = c(b1,b2)

	#### Establish color palettes
	#### NB Change white to a whitish shade
	
	c1 = colorRampPalette(c('purple','blue','green'))
	c2 = colorRampPalette(c('yellow','orange','red'))

	### Get a list of colors based on the length of breaks from the color palettes

	cd1 = c1(length(b1))
	cd2 = c2(length(b2))
	allcodes = c(cd1,'#FFFFFF',cd2) # This a list of colors one element longer than all breaks, with white lining up with zero
	
	### Open the .png device driver
	
	png(paste(plot.dir,spp,'-Diff.png',sep=''),units='cm',height=10,width=10,res=1000)

	# Set plotting parameters
	
	par(oma=c(0,.5,.5,0), bg='#FFFFFF00')
	
	# Plot DEM

	image(dmap, col=allcodes, bty='n', axes=F, xlab='', ylab='')
	
	### Use a for loop to plot a horizontal legend gradient

	### Establish x-y's for the corners

	pnts = data.frame(x=c(144.7,144.7,144.9,144.9),y=c(-17,-19,-19,-17))

	### Establish break points for gradient plotting

	yvals = seq(min(pnts[, 2]), max(pnts[, 2]), (max(pnts[, 2]) - min(pnts[, 2]))/(length(allcodes)+1))
		 
	for (inc in 1:length(allcodes))

		{
		
		 polygon(x=pnts[,1],y=c(yvals[inc],yvals[inc+1],yvals[inc+1],yvals[inc]),col=allcodes[inc],border=NA)
		
		}

	### Plot a hollow polygon with a thin black border around the DEM scale

	polygon(x=pnts$x,y=c(max(yvals),min(yvals),min(yvals),max(yvals)),col='#FFFFFF00',lwd=.5)
	
	#### Text labels
	
	text(144.8,-16.9,paste(round(max(lims,na.rm=T),1)),cex=.6)
	text(144.8,-19.1,paste(round(min(lims,na.rm=T),1)),cex=.6)
	
	dev.off()
	
	cat('\n',spp,'\n')
	
	
	}
	
# Done
	
		
### Begin plotting, first loop by stat type, then species plotting summaries of ClassStat Data
### Make a colour palette first

#### Establish color palette
	
cols = colorRampPalette(c('black','blue','purple','red'))

### Start looping through stats

for (stat in c(4:13,16:18,21:22,25:26,29:32,35:41))

	{
	
	### subset data to a single statistic
	
	tdata = outdata[c(1:3,stat)]
	
	### Open .png device driver
	
	png(paste(plot.dir,names(outdata)[stat],'.png',sep=''),width=18,height=18,units='cm',res=1000)
	
	### Establish plotting limits
	
	lims = range(tdata[,4])
	
	### Configure the plot space
	
	plot(tdata[,4],tdata[,4],xlim=lims,ylim=lims,main=names(outdata)[stat],xlab='BIOCLIM',ylab='accuCLIM', type='n')
	
	### Track changes
	
	i=1
	
	### Add a legend
	
	legend('bottomright',legend=(unique(tdata[,1])), col=cols(7), text.col = cols(7), pch=c(1:7), bty='n')
	
	### Add a 1:1 line
	
	abline(a=0,b=1,col='black',lty=2,lwd=1.2)
	
	### Loop through species
	
	for (spp in unique(tdata[,1]))
	
		{
		
		### Subset to a single species
		
		ttdata = tdata[which(tdata[,1]==spp),]
		
		### Plot the points
		
		points(ttdata[1,4],ttdata[2,4],col=cols(7)[i], cex=1.6,pch=i,type='p')
		
		### Track changes
		
		i=i+1
		
		### Report progress
		
		cat('\n',stat,'\n')
		
		}
		
	dev.off()
	
	}
	
# Done

### Now plot up raw MaxEnt ASCII's against one another

for (spp in species)

	{
	
	### Read in the ASCIIs
		
	bcmap = read.asc.gz(paste(asc.dir,spp,'/BIOCLIM/output/',spp,'_BIOCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	acmap = read.asc.gz(paste(asc.dir,spp,'/microCLIM/output/',spp,'_microCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	
	### Remove NA's from maps
	
	bcmap = bcmap[which(is.na(bcmap)==F)]
	acmap = acmap[which(is.na(acmap)==F)]
		
	### Establish plotting limits
	
	lims = range(c(acmap,bcmap))
	
	### Plot the data
		
	png(paste(plot.dir,spp,'RawMaxEntComparison.png',sep=''),height=18,width=18, units='cm',res=1000)
	
	### Plot
	
	plot(bcmap,acmap, xlab='BIOCLIM ES', ylab = 'accuCLIM ES',main=spp, col='#FF000075', pch=1, xlim=c(0,1), ylim=c(0,1))
	
	### 1:1 line
	
	abline(a=0,b=1,col='black',lty=2,lwd=.9)
	
	### Turn off device
	
	dev.off()
	
	### Display progress
	
	cat('\n',spp,'\n')
	
	}
	

############### Looking at distribution of accuCLIM ES values within bins of BIOCLIM ES values

for (spp in species)

	{
	
	### Read in the ASCIIs
		
	bcmap = read.asc.gz(paste(asc.dir,spp,'/BIOCLIM/output/',spp,'_BIOCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	acmap = read.asc.gz(paste(asc.dir,spp,'/microCLIM/output/',spp,'_microCLIM_Realized_No_Small_Patches.asc.gz',sep=''))
	
	### Establish positions
	
	base.pos = as.data.frame(which(is.finite(bcmap), arr.ind=T))
	
	base.pos$lat = getXYcoords(bcmap)$x[base.pos$row]
	base.pos$long = getXYcoords(bcmap)$y[base.pos$col]
	
	### Extract data from ASCII's
	
	base.pos$ac = extract.data(cbind(base.pos$lat,base.pos$long),acmap)
	base.pos$bc = extract.data(cbind(base.pos$lat,base.pos$long),bcmap)
	
	### Set some plotting parameters and start the .png device driver
	
	png(paste(plot.dir,spp,'-accuCLIM_Frequency_Panels.png',sep=''),units='cm',width=20,height=25,res=1200)
	
	par(mfrow=c(4,5))

	### Start looping by percentile to establish max y-limit for plotting
	
	for (bin in seq(0,.95,.05))
	
		{
		
		### Select data from base.pos within a single bin of BIOCLIM ES values
	
		t.pos = base.pos[which(base.pos$bc>bin & base.pos$bc<=bin+.05),]
		
		if(nrow(t.pos)>1)
		
		{
		
		### Create a density object of accuCLIM ES values to plot
		
		toplot = density(t.pos$ac[which(is.finite(t.pos$ac)==T)])
		
		### Count up points below and above expected values
		
		below = t.pos[which(t.pos$ac>=bin & t.pos$ac<((bin+bin+.05)/2)),]
		above = t.pos[which(t.pos$ac>=((bin+bin+.05)/2) & t.pos$ac<(bin+.05)),]
		
		### Plot density object
		
		plot(toplot,xlab='accuCLIM ES',ylab='ECDF',xlim=c(0,1),main=paste(spp,' BC ES ',bin,'-',bin+.05,sep=''), lwd=1.2,col='lightblue',cex.main=.7)
		
		### Put on a line representing the expected value
		
		lines(x=rep((bin+bin+.05)/2,2),y=c(0,max(toplot$y)),lty=2,lwd=1.5,col='lightblue',type='l')
		
		### Put on a legend
		
		legend('bottomright',legend=c(paste(nrow(above))),text.col='red',bty='n',cex=.8)
		legend('topleft',legend=c(paste(nrow(below))),text.col = 'darkblue',bty='n',cex=.8)
		
				
		}
		
		}
		
	### Turn off device after writing all 20 plots to the same .png file
		
	dev.off()
		
	### Report progress
		
	cat('\n',spp,'\n')
		
	}
	
### Done

################## Create maps of distributions

for (spp in species)

	{
	
	### Loop through different models within species, extracting and summarizing data
	
	for (mod in c('BIOCLIM','microCLIM'))
	
		{
		
		### Read in the ASCII
		
		t.asc = read.asc.gz(paste(asc.dir,spp,'/',mod,'/output/',spp,'_',mod,'_Realized_No_Small_Patches.asc.gz',sep=''))
		
		### Create a color palette
		
		c1 = colorRampPalette(c('white','yellow','orange','red'))
	
		### Open the .png device driver
	
		png(paste(plot.dir,spp,'-',mod,'-Distribution.png',sep=''),units='cm',height=12.5,width=10,res=1000)

		# Set plotting parameters
	
		par(bg='lightgrey')
	
		# Plot DEM

		image(t.asc, col=c1(100), bty='n', axes=F, xlab='', ylab='')
		
		dev.off()
		
		cat('/n',spp,'-',mod,'\n')
		
		}
		
	}
	
### Done
	
		
		
		
	
		
	
		

		
	
	
	
	
	
	
	
